MODULE INIT
  USE GLOBVAR
  !USE MPI
  !USE INTEGRATION
  IMPLICIT NONE
  
CONTAINS
  
  SUBROUTINE INITIALIZE()
    IMPLICIT NONE
    ! data fill in 
    INTEGER, PARAMETER :: N_vars=50     ! Number of columns in dataset
    REAL(8) :: dataset(N_ind*N_test,N_vars),pi
    INTEGER :: i,j,k,t,m,l
    
    ! data moment
    !REAL(8) :: p(N_testl,N_test)
    !INTEGER :: mk(N_testl,N_test),nk(N_testl,N_test)




    ! Now read in the dataset and fill in the data matrices
    IF (myid==0) OPEN(176,file="learning_task_fortran_0617.raw")
    DO i=1,N_ind*N_test
       IF (myid==0) READ(176,fmt=*) dataset(i,:)
    END DO
       
    DO m=1,N_ind*N_test
       i=INT((m-1)/N_test)+1
       l=m-((i-1)*N_test)
       
       ! fill in language tasks
       DO j=1,N_testl
           IF (INT(dataset(m,15))==j) THEN
             l_test(i,j)=dataset(m,6)        ! task performance
             IF (l_test(i,j)>-1) m_l(i,j)=1  ! Whether observe the performance
             IF (l_test(i,j)<-1) m_l(i,j)=-1
             x_testl(i,j,1:5)=dataset(m,7:11)
             x_testl(i,j,6)  =1.0d0
             indl(i,j)=dataset(m,31)
             repl(i,j)=dataset(m,35)
             levell(j)=dataset(m,23)
             repeatl(j)=dataset(m,27)
             monthl(j)=dataset(m,4)
             timel(j)=dataset(m,5)
           END IF
           IF (INT(dataset(m,19))==j) THEN
               firstl(i,1)=j
               firstl(i,2)=dataset(m,3)
           END IF
       END DO
       
       ! fill in cognitive tasks
       DO j=1,N_testc
           IF (INT(dataset(m,16))==j) THEN
             c_test(i,j)=dataset(m,6)        ! task performance
             IF (c_test(i,j)>-1) m_c(i,j)=1  ! Whether observe the performance
             IF (c_test(i,j)<-1) m_c(i,j)=-1
             x_testc(i,j,1:5)=dataset(m,7:11)
             x_testc(i,j,6)  =1.0d0
             indc(i,j)=dataset(m,32)
             repc(i,j)=dataset(m,36)
             levelc(j)=dataset(m,24)
             repeatc(j)=dataset(m,28)
           END IF

           IF (INT(dataset(m,20))==j) THEN
               firstc(i,1)=j
               firstc(i,2)=dataset(m,3)
           END IF
           
       END DO
       
       ! fill in fine motor tasks
       DO j=1,N_testf
           IF (INT(dataset(m,17))==j) THEN
             f_test(i,j)=dataset(m,6)        ! task performance
             IF (f_test(i,j)>-1) m_f(i,j)=1  ! Whether observe the performance
             IF (f_test(i,j)<-1) m_f(i,j)=-1
             x_testf(i,j,1:5)=dataset(m,7:11)
             x_testf(i,j,6)  =1.0d0
             indf(i,j)=dataset(m,33)
             repf(i,j)=dataset(m,37)
             levelf(j)=dataset(m,25)
             repeatf(j)=dataset(m,29)
           END IF
           
           IF (INT(dataset(m,21))==j) THEN
               firstf(i,1)=j
               firstf(i,2)=dataset(m,3)
           END IF
       END DO
       
        ! fill in gross motor tasks
       DO j=1,N_testg
           IF (INT(dataset(m,18))==j) THEN
             g_test(i,j)=dataset(m,6)        ! task performance
             IF (g_test(i,j)>-1) m_g(i,j)=1  ! Whether observe the performance
             IF (g_test(i,j)<-1) m_g(i,j)=-1
             x_testg(i,j,1:5)=dataset(m,7:11)
             x_testg(i,j,6)  =1.0d0
             indg(i,j)=dataset(m,34)
             repg(i,j)=dataset(m,38)
             levelg(j)=dataset(m,26)
             repeatg(j)=dataset(m,30)
           END IF
           
            IF (INT(dataset(m,22))==j) THEN
               firstg(i,1)=j
               firstg(i,2)=dataset(m,3)
           END IF
       END DO
       
       IF (l==1) THEN
           L0(i,1:4)=dataset(m,11:14)
           L0(i,5)  =dataset(m,43)
           L0(i,6)  =1.0d0
           intermean(i)=dataset(m,50)
                      id(i)=dataset(m,1)
       END IF
    END DO
    
    IF (myid==0) THEN 
       CLOSE(176)
       WRITE(*,*) "dataset read"
    END IF
    !!!!!!!!!!!!!!!!!! Calculate Data Moment!!!!!!!!!!!!!!!!!!!!!!!!!
    !
    !DO i=1,N_testl*N_test
    !    k=INT((i-1)/N_test)+1
    !    j=i-(k-1)*N_test
    !    !write(*,*) i,k,j
    !    moment(i)=DBLE(SUM(y_test(:,k,j))/DBLE(N_IND))
    !END DO
    !
    !!!!! Need More Moments
    !
    !!!! Calculate conditional probability 
    !    mk=0
    !    nk=0
    !    DO j=1,N_testl
    !        DO k=2,N_test
    !            DO i=1,N_ind
    !                IF (y_test(i,j,k-1)==1) THEN
    !                    mk(j,k-1)=mk(j,k-1)+1
    !                    IF (y_test(i,j,k)==1) THEN
    !                       nk(j,k-1)=nk(j,k-1)+1
    !                    END IF
    !                END IF
    !            END DO
    !        END DO
    !    END DO
    !    DO j=1,N_testl
    !        DO k=1,N_test-1
    !            p(j,k)=DBLE(DBLE(nk(j,k))/DBLE(mk(j,k)))
    !        END DO
    !    END DO
    !    
    !DO i=N_testl*N_test+1,N_testl*(2*N_test-1)
    !    k=INT((i-N_testl*N_test-1)/(N_test-1))+1
    !    j=i-N_testl*N_test-(k-1)*(N_test-1)
    !    !write(*,*) i,k,j
    !    moment(i)=p(k,j)
    !END DO         
    !    mk=0
    !    nk=0
    !    p=0.0d0
    !    DO j=1,N_testl
    !        DO k=3,N_test
    !            DO i=1,N_ind
    !                IF (y_test(i,j,k-2)==1) THEN
    !                    mk(j,k-2)=mk(j,k-2)+1
    !                    IF (y_test(i,j,k)==1) THEN
    !                       nk(j,k-2)=nk(j,k-2)+1
    !                    END IF
    !                END IF
    !            END DO
    !        END DO
    !    END DO
    !    DO j=1,N_testl
    !        DO k=1,N_test-2
    !            p(j,k)=DBLE(DBLE(nk(j,k))/DBLE(mk(j,k)))
    !        END DO
    !    END DO
    !DO i=N_testl*(2*N_test-1)+1,N_testl*(3*N_test-3)
    !    k=INT((i-N_testl*(2*N_test-1)-1)/(N_test-2))+1
    !    j=i-N_testl*(2*N_test-1)-(k-1)*(N_test-2)
    !    write(*,*) i,k,j
    !    moment(i)=p(k,j)
    !END DO     
    !
    !    mk=0
    !    nk=0
    !    p=0.0d0
    !    DO j=2,N_testl
    !        DO k=1,N_test
    !            DO i=1,N_ind
    !                IF (y_test(i,j-1,k)==1) THEN
    !                    mk(j-1,k)=mk(j-1,k)+1
    !                    IF (y_test(i,j,k)==1) THEN
    !                       nk(j-1,k)=nk(j-1,k)+1
    !                    END IF
    !                END IF
    !            END DO
    !        END DO
    !    END DO
    !    DO j=1,N_testl-1
    !        DO k=1,N_test
    !            p(j,k)=DBLE(DBLE(nk(j,k))/DBLE(mk(j,k)))
    !        END DO
    !    END DO
    !DO i=N_testl*(3*N_test-3)+1,N_testl*(3*N_test-3)+(N_testl-1)*N_test
    !    k=INT((i-N_testl*(3*N_test-3)-1)/(N_test))+1
    !    j=i-N_testl*(3*N_test-3)-(k-1)*(N_test)
    !    !write(*,*) i,k,j
    !    moment(i)=p(k,j)
    !END DO     
    !    mk=0
    !    nk=0
    !    p=0.0d0
    !    DO j=1,N_testl
    !        DO k=2,N_test
    !            DO i=1,N_ind
    !                IF (y_test(i,j,k-1)==0) THEN
    !                    mk(j,k-1)=mk(j,k-1)+1
    !                    IF (y_test(i,j,k)==0) THEN
    !                       nk(j,k-1)=nk(j,k-1)+1
    !                    END IF
    !                END IF
    !            END DO
    !        END DO
    !    END DO
    !    DO j=1,N_testl
    !        DO k=1,N_test-1
    !            p(j,k)=DBLE(DBLE(nk(j,k))/DBLE(mk(j,k)))
    !        END DO
    !    END DO   
    !    
    !DO i=N_testl*(3*N_test-3)+(N_testl-1)*N_test+1,N_testl*(3*N_test-3)+(N_testl-1)*N_test+N_testl*(N_test-1)
    !    k=INT((i-N_testl*(3*N_test-3)-(N_testl-1)*N_test-1)/(N_test-1))+1
    !    j=i-N_testl*(3*N_test-3)-(N_testl-1)*N_test-(k-1)*(N_test-1)
    !    !write(*,*) i,k,j
    !    moment(i)=p(k,j)
    !END DO 
    !
    !
    !    mk=0
    !    nk=0
    !    p=0.0d0
    !    DO j=1,N_testl
    !        DO k=3,N_test
    !            DO i=1,N_ind
    !                IF (y_test(i,j,k-2)==0) THEN
    !                    mk(j,k-2)=mk(j,k-2)+1
    !                    IF (y_test(i,j,k)==0) THEN
    !                       nk(j,k-2)=nk(j,k-2)+1
    !                    END IF
    !                END IF
    !            END DO
    !        END DO
    !    END DO
    !    DO j=1,N_testl
    !        DO k=1,N_test-2
    !            p(j,k)=DBLE(DBLE(nk(j,k))/DBLE(mk(j,k)))
    !        END DO
    !    END DO   
    !    
    !DO i=N_testl*(3*N_test-3)+(N_testl-1)*N_test+N_testl*(N_test-1)+1,N_testl*(3*N_test-3)+(N_testl-1)*N_test+N_testl*(2*N_test-3)
    !    k=INT((i-(N_testl*(3*N_test-3)+(N_testl-1)*N_test+N_testl*(N_test-1))-1)/(N_test-2))+1
    !    j=i-(N_testl*(3*N_test-3)+(N_testl-1)*N_test+N_testl*(N_test-1))-(k-1)*(N_test-2)
    !    !write(*,*) i,k,j
    !    moment(i)=p(k,j)
    !END DO 
    !
    !    mk=0
    !    nk=0
    !    p=0.0d0
    !    DO j=2,N_testl
    !        DO k=1,N_test
    !            DO i=1,N_ind
    !                IF (y_test(i,j-1,k)==0) THEN
    !                    mk(j-1,k)=mk(j-1,k)+1
    !                    IF (y_test(i,j,k)==0) THEN
    !                       nk(j-1,k)=nk(j-1,k)+1
    !                    END IF
    !                END IF
    !            END DO
    !        END DO
    !    END DO
    !    DO j=1,N_testl-1
    !        DO k=1,N_test
    !            p(j,k)=DBLE(DBLE(nk(j,k))/DBLE(mk(j,k)))
    !        END DO
    !    END DO   
    !    
    !DO i=N_testl*(3*N_test-3)+(N_testl-1)*N_test+N_testl*(2*N_test-3)+1,N_testl*(3*N_test-3)+2*(N_testl-1)*N_test+N_testl*(2*N_test-3)
    !    k=INT((i-(N_testl*(3*N_test-3)+(N_testl-1)*N_test+N_testl*(2*N_test-3))-1)/(N_test))+1
    !    j=i-(N_testl*(3*N_test-3)+(N_testl-1)*N_test+N_testl*(2*N_test-3))-(k-1)*(N_test)
    !    !write(*,*) i,k,j
    !    moment(i)=p(k,j)
    !END DO
    
                !OPEN(6424,file="ltest.out")
                !  DO i = 1,N_ind
                !    WRITE(6424,'(200I4)') (l_test(i,j),j=1,N_testl)
                !  END DO
                !CLOSE(6424)
                !OPEN(6424,file="indl.out")
                !  DO i = 1,N_ind
                !    WRITE(6424,'(200I4)') (indl(i,j),j=1,N_testl)
                !  END DO
                !CLOSE(6424)
                !OPEN(6424,file="repl.out")
                !  DO i = 1,N_ind
                !    WRITE(6424,'(200I4)') (repl(i,j),j=1,N_testl)
                !  END DO
                !CLOSE(6424)
                !OPEN(6424,file="data.out")
                !  DO i = 1,N_ind*N_test
                !    WRITE(6424,'(200F16.8)') (dataset(i,j),j=1,N_vars)
                !  END DO
                !CLOSE(6424)
    
    


    ! For integration
    !CALL Gauss_Hermite(xgh,wgh)
    !pi=3.14159265359d0
    !wgh=(pi**(-0.50d0))*wgh

  END SUBROUTINE INITIALIZE
  
END MODULE INIT
